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We study the ground state properties of an ABA-stacked trilayer graphene. The low energy band 
structure can be described by a combination of both a linear and a quadratic particle-hole sym- 
metric dispersions, reminiscent of monolayer- and bilayer-graphene, respectively. The multi-band 
structure offers more channels for instability towards ferromagnetism when the Coulomb interaction 
is taken into account. Indeed, if one associates a pseudo-spin 1/2 degree of freedom to the bands 
(parabolic/linear), it is possible to realize also a band-ferromagnetic state, where there is a shift 
in the energy bands, since they fill up differently. By using a variational procedure, we compute 
the exchange energies for all possible variational ground states and identify the parameter space for 
the occurrence of spin- and band-ferromagnetic instabilities as a function of doping and interaction 
strength. 
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I. INTRODUCTION 

The successful isolation of a one atom thick carbon 
layer, graphene, has attracted enormous interest in the 
field of condensed matter.EEl One intriguing aspect of 
the problem is that upon coupling a finite number of 
graphene layers, novel and unexpected properties emerge. 
Compared to the strong sp^ bonding between carbon 
atoms within the graphene sheet, the weak van der Waals 
force between the layers allows for the formation of dif- 
ferent hybridized iV-layered configurations. The resulting 
system is then different from both its 2D (graphene) and 
its 3D (graphite) counterparts, and depends strongly on 
the number of layers and on how the stacking is realized. 
The investigation of multi-layer graphene may open new 
avenues in the understanding of graphene's electronic 
properties and in the field of device engineering 

Many of the unique electronic properties of monolayer 
graphene, as opposed to the more conventional GaAs 
2D electron gas, originate from the geometry of the 
honeycomb lattice. These include the peculiar gapless 
Dirac-cone dispersion,'^! the unconventional integer quan- 
tum Hall effect,'^ and Klein tunneling,^ to name a few. 
On the other hand, multi-layer graphene exhibits differ- 
ent but equally interesting features. While the particle- 
hole symmetry is generally preserved in the band struc- 
ture obtained from the minimal tight-binding descrip- 
tion, the number of conical points and the low-energy 
dispersion both depend sensitively on the stacking con- 
figuration of the iV-layered structure. For example, in 
the so-called Bernal stacking of a bilayer graphene, the 
conduction and valence bands touch at the same two 
points in the Brillouin zone as they do in monolayer 
graphene, but disperse quadratically instead of linearly. 
This feature has attracted much interest because it allows 
for strong electron correlations to take place.l^ Very re- 
cently, broken-symmetry states have been observed due 
to interaction effects in suspended bilayer graphene .'^''^ 
Although a complete characterization of their proper- 



ties is still lacking, there are some interesting theoreti- 
cal proposals for the observed states: a many-body exci- 
tonic stat^^ or a n anom alous spin-Hall state with time- 
reversal symmetry l^^'^'^' For another example, the relative 
twist angle in a bilayer graphene can lead to a highly com- 
plex Moire band structure, which requires a description 
beyond the standard Bloch's band picture. In fact, at a 
particular twisting angle, the van Hove singularity of the 
usual graphene band structure can become observable at 
a relatively low energy of a few meV.^ Since high quality 
samples of A^-layered graphene are now becoming acces- 
sible experimentally, their anticipated new properties are 
just about to be unraveled. 

In trilayer graphene, the transport properties are also 
different, depending on the stacking order: at the Dirac 
point, the ABA-stacked trilayer (Fig. [T]) is a semimetal, 
whereas the ABC one is a semiconductor, with an intrin- 
sic band gap.l^ The electronic band structure in ABC- 
stacked trilayer graphene was determined using an effec- 
tive mass approximatiorP^ and using an ab-initio den- 
sity functional theory.l^ On the other hand, for ABA- 
stacked the band structure was calculated in the pres- 
ence of external gates using a self-consistent Hartree 
approximation.E^In the absence of a gate, the low-energy 
spectrum consists of superimposed linear and quadratic 
bands, which touch at k = 0. In the presence of a mag- 
netic field, the plateau structure in the Hall conductivity 
is also determined by the stacking order. Very recently, 
the integer quantum Hall effect was e xperimentally ob- 
served in an ABC-stacked sample.^^It was shown that 
the effect is similar to the one observed in monolayer 
graphene,'^ except for the first plateau at filling factor 
V = 2^ which was not observed in the trilayer sample. 
Indeed, this plateau is governed by the chirality of the 
quasi-particles, which is 1, 2, and 3 for monolayer, bi- 
layer, and trilayer graphene, respectively. The corre- 
sponding Berry phases are thus tt, 27r, and 37r, respec- 
tively. With regard to the ABA stacking, the problem 
of low mobility has been recently overcome, by grow- 
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Figure 1. (Color online) ABA stacked trilayer graphene with 
the various hopping parameters. 



ing the sample on a high-quality hexagonal boron nitride 
substrate, which reduces the carrier scattering. The pe- 
culiar crossing of the Landau levels due to the massive 
and massless sub-bands has allowed for a direct deter- 
mination of the Slonczweski-Weiss-McClure model pa- 
rameters used to describe the electronic structure of the 
material!^^ 

We focus here on the ground state properties of tri- 
layer graphene in the ABA-stacking configuration in the 
presence of interactions. The ground state of A^-layer 
undoped graphene is usually assumed to be the state 
in which the energy bands are filled up to the Dirac 
point. However, the energy bands are spin degenerate 
and the formation of pockets of opposite sign in the two 
spin degenerate bands leads to a gain in exchange en- 
ergy. This gain in exchange energy is accompanied by a 
cost in kinetic energy. In monolayer graphene, the cost 
in kinetic energy is large enough to prevent any ferro- 
magnetic instabilitiesj231 only if the interaction would be 
tuned to unphysical values one would observe the spon- 
taneous generation of spin up and spin down pockets. 
In bilayer graphene the situation is different. The lead- 
ing order term in the exchange energy is one order lower 
in the pocket size than the kinetic energy is. Therefore, 
the exchange interaction dominates and pockets will form 
with a size, in fc-space, of order Q sa 0.05t±, where tj_ is 



the interlayer hopping energy in dimensionless units and 
Q is measured in units of some cut off."^ Hence, bilayer 
graphene has a small ferromagnetic instability. The coex- 
istence of a parabolic and a linear bands in ABA-trilayer 
graphene opens the way to investigate, next to ordinary 
ferromagnetic instabilities (Fig. [2^), also the 'band ferro- 
magnetism' phenomenon. With band ferromagnetism we 
mean that the two bands (linear and parabolic) become 
shifted with respect to each other (the crossing point of 
the linear and parabolic conduction and valence bands no 
longer overlap) , or alternatively, that the bands fill up to 
different Fermi energies (see Fig.[2]D). In the following, we 
will generalize the approach used in Refs.|23|and[2S]to in- 
vestigate ferromagnetic instabilities in trilayer graphene. 
We will show that spin- and band-ferromagnetism may 
occur both separately and simultaneously. The paper is 
organized as follows: In section II we introduce the model 
that we use in section III to compute (band) ferromag- 
netic instabilities for both undoped and doped trilayer 
graphene. Our conclusions are presented in section IV. 
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Figure 2. (Color online) Sketch of (a) the spin- ferromagnetic 
state in an undoped trilayer and (b) the band-ferromagnetic 
state in a doped trilayer. 



II. THE MODEL 

In this paper, we use a tight-binding approximation 
to model trilayer graphene and perform an expansion 
around the K point. The low-energy Hamiltonian around 
the K point is given by 



p,cn 
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and the sum is over all relevant quantum numbers. Here, 
p a i^ipa) creatcs a particle with momentum p and 
spin a at the A [B) sublattice in the i-th layer {i = 



1,2,3), t_L ~ 0.35 eV is the interlayer hopping energy, 
vf — {3/2)at denotes the Fermi velocity, with a — 0.142 
nm the lattice spacing and i « 3 eV the nearest neighbor 
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hopping energy, p is the norm of the momentum vector 
P = {Px,Py) and = aictan {py / px) ■ Note that if one 
would have expanded around the K' point, we would 
have found a Hamiltonian which is the complex conju- 
gate of Eq. ([!]). Since we neglect intervalley interactions, 
we do not need to take this into account and we simply 
multiply our results by a factor two. 

We perform a change of basis, 4" with 
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to bring the Hamiltonian into the form H 
E*J,,,H(p)^p,,, where 
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Thus, the trilayer can be described as a combination of a 
monolayer and a bilayer with a modified interlayer hop- 
ping energy. Note that in the new basis, the basis vectors 
that are associated with the monolayer part are odd un- 
der reflection with respect to the middle plane, while the 
ones that describe the bilayer are even under this trans- 
formation. The hopping parameters 72 and 75 from the 
Slonczewski-Weiss-McClure (SWM)-model, or a voltage 
difference between the top and bottom layer break this 
reflection symmetry and couple the blocks in the trilayer 
Hamiltonian.'^ We will neglect those terms here. 

Since the Hamiltonian has a block form and we know 
how to diagonalize the different blocks, it is now a trivial 
task to bring it into a diagonal form. Using the results 
from Refs. [24l and l25l we find that H{p) can be diago- 
nalized as follows: 



P(p) = W\p)Hip)W{p) 
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Figure 3. (Color online) The energy spectrum of tri- 
layer graphene. The numbering of the bands is such that 



where V"(p) and M(p) are the matrices that diagonalize 
the monolayer and bilayer Hamiltonian respectively. 
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In the last matrix, <y3(p) is defined by the relation 
tan[2(/7(p)] = VFy/^plt±. This result differs by a fac- 
tor ^/2 from Ref. because of the modified interlayer 
hopping parameter in T-Lu- The energy bands are given 
by the nonzero entries of the matrix I'(p), 



X'(p) = diag{ - vfP, vpp, 
[-t±+ap)]/V2, [tj_ 



where ^ (p) = + 2vpp'^. 

The next step is to implement the Coulomb interac- 
tion in the model. Since we consider only weakly doped 
trilayers in this paper, the Coulomb interaction is only 
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slightly screened and therefore long ranged, 

Hi^lJ d2xd2y{^°(x-y)[pi(x)pi(y)+p2(x)p2(y) 

+ p3 (x)p3 (y)] + (x - y) [pi (x)p2 (y) 

+ P2(x)/5l(y) + P2(x)p3(y) + P3(x)p2(y)] 

+ F^^D (x - y) [p, (x)p3 (y) + p3 (x)pi (y)] } , (2) 

where p,(x) = E<t (^U W"*^'^ W + W^'^'^ W) 
the density of electrons in the i-th layer and the interac- 
tion potentials for the in-plane (D), the nearest-neighbor 
planes (ND) and the next-nearest-neighbor planes (2ND) 
are given by 
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Here, e is the electron charge, e the dielectric constant of 
the substrate (of air in the case of suspended graphene), 
and d the interlayer distance {d « .32 nm). The form 
of (x — y) can be understood by recalling that x is 
a 2 dimensional vector. We Fourier transform Eq. ([2| 
and express it in terms of symmetric and anti-symmetric 
combinations of layer densities. 



[Pa{ci)Va{q)pai-q) 



2yl^ ^ 

q a—± 

+ PQ(q)Va(q)pa(-q) + Pa(q)t4(q)pa(~q)], (3) 

where the prime on the sum indicates that we omit the 
q = term, since it is canceled by the neutralizing back- 
ground (Jellium model), A is the area of the unit cell, 
and the different quantities are defined by 
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We want to write this interaction term in the number 
operators of the energy bands instead of the number op- 
erators of the layers. We know how to diagonalize the 
kinetic term and therefore ^p^a = ■^(p)^^'p.cr are the op- 
erators that annihilate particles in the different energy 



bands. As a result, we obtain, (<I>J, ,^$p_cr)^ = "^.^(p): 
the number operator of the j-th energy band, where we 
have to number the bands as in Fig. [3] It is convenient 
to rewrite the density operators in the diagonal basis, 

p±(q) = 5I*I.+q^*(p + q'P)*p' (9) 

p 

p±(q) = 5I*I>+q^*(p + q'P)*p' (10) 

p 

p±(q) = ^*p+qX*(p + q,p)$p, (11) 



where 

x'^lp + q^p) 
x'^(p + q,p) 
x'^(p + q,p) 



— Z^+qdiag(l,l,±l,±l,0,0)Zp, (12) 



:ZT+qdiag(0,0,±l,±l,l,l)Zp, (13) 
:Z^+qdiag(l,l,0,0,±l,±l)Zp. (14) 



Inserting equations ([4|-(14) into the interaction Hamil- 
tonian ([3| yields the interaction term that we use for 
our calculations. We are only interested in the exchange 
energy, which is given by 
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Xrj(p',P)Xjj(P,P')^a(p' - P)nz.a,a{P'hj.<yAP) 

Xrj(p',P)Xj»(P,P')K»(p' - P)n^.aAp')nj.aAP) 

(15) 

In the sum a takes the values ±; i and j label compo- 
nents, hence run from 1 to 6; cr sums over spin, and a 
over the valley index. We neglected the valley index so 
far since in our case it only gives rise to an extra fac- 
tor two, as we choose the same pocket structure for both 
valleys in our studies. 



III. FERROMAGNETIC INSTABILITIES 

Undoped case 

For undoped trilayer graphene, the noninteracting 
ground state is the configuration in which the three va- 
lence bands are completely filled and the conduction 
bands are completely empty. If an electron or hole pocket 
forms in one of the bands, this costs kinetic energy. 
This cost is given by the absolute value of the integral 
/p^^'^' dE p{E)E, where Q is the pocket size and p{E) the 
density of states. Since for the linear band, p{E) ~ E 
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parameters. For zero doping one has the constraint: 




Figure 4. (Color online) The energy 

AE{Qiu,Qid,Qpu,Qpd) per unit cell (Eq 
undoped trilayer, where we have chosen Q;„ 
Because of particle number conservation Qp„ — —Qpd ^ 
AE is measured in units of hvpA, Qi and Qp are 
measured in units of A. 



difference 
18|) for the 

— Qld ~ ^ 



both 



and E{ff) ^ Q, one finds that Ai?kin,; ^ ^ while for 
the parabolic band, p{E) ~ but E{ff) ~ Q^, hence 
Ai?kin,p ^ Q"^- In fact, the changes in kinetic energy for a 
linear band with pockets of size Qi and a parabolic band 
with pockets of size Qp are 



A£;ki„,i(Qi) 



A 



Ai?ki„,p(Qp)-^^^|Qpr 

Since Qi < 1, for i = ^/p, the creation of linear pockets 
costs more kinetic energy than the creation of parabolic 
ones. Trilayer graphene has four energy bands close to 
the K point, hence there are four different pocket param- 
eters: Qiu, Qld, Qpu, and Qpd, where l/p stands for linear 
and parabolic bands and u/d for up and down spins. We 
are assuming long range interactions and are neglecting 
the short range part, hence there is no intervalley scat- 
tering. We also assume particle number conservation, 
thus Qpd is not independent from the other variational 
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Figure 5. (Color online) The energy difference per unit cell 
AE{Qiu = 0, Qld = 0, Qpu = Qp, Qpd = —Qp)- This is a cross 
section along the Qi — axis of Fig. |4] AE is measured in 
units of hvpA and Qp is measured in units of A. 
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(16) 



where Sia- = 4-1 for electron-like pockets and Si^ = — 1 
for hole-like pockets. 

One can now vary the pocket parameters and calculate 
whether the energy is minimized for nonzero pocket sizes 
(at zero temperature). Our formalism is build up in such 
a way that the pocket parameters can take both positive 
and negative values. A positive Q corresponds to an elec- 
tron pocket. Hence, the corresponding conduction band 
(linear/parabolic, up/down) is filled up to momentum Q. 
A negative Q corresponds to hole pockets, i.e. the cor- 
responding valence bands are depleted up to momentum 
\Q\. This method allows us to obtain the exchange inte- 
grals for all possible pocket configurations at once. Using 
this formalism, we find that the bands fill up according 
to (see Fig. |3] for the numbering of the bands) , 
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where 8 is the Heaviside step function. Note that one 
cannot have both electron and hole pockets in the same 
band at the same time, because if, for example, Qiu > 0, 
then Q{—Qiu — p) = 0. Hence, in this case, the lin- 
ear spin up valence band is completely filled (band 1 in 
Fig. [3]), while the linear spin up conduction band (band 2 
in Fig. [3]) is filled up to momentum corresponding 
to an electron pocket of size |(5;u|. 

The integrals that we have to compute have the same 
structure as the ones in Ref. [551 The expansion in the 
pocket parameters is highly nontrivial and very lengthy. 
Since there are three variational parameters, we have per- 
formed the integrals numerically. The expression for the 
integrals (Eq. [Tsl ) has many terms and it is not enlight- 
ening to write all of them out. 

From this point on, we work in dimensionless units by 
measuring momenta in units of a cutoff A, which is esti- 
mated using a Debye approximation, in which the num- 
ber of states is conserved in the Brillouin zone: A^ = 
2Tr/A. We measure energies in units of hvpA{AA'^) = 
hvp^- This dimensionless energy corresponds with the 
energy per unit cell in units of hvpA. Let us also intro- 
duce a dimensionless interaction strength g — j (ehvp)- 
Furthermore, we set A, h and t equal to unity. Note that 
the spin-up and spin-down terms decouple. This allows 



6 



us to calculate 

AE{Qi,Qp) = AEkin{Qi,Qp) + AE,^{Qi,Qp) 

+ AE^^,i{Qi) + AE^^^p{Qp) 

.mixed 

{Qi,Qp) (17) 

on a discrete Ni x Np lattice, where we have chosen the 
values of the pocket parameters such that their squares lie 
on an equally spaced grid for reasons which will become 
clear later. After calculating these data points, one can 
compute 

Qpd) — AE{Qiu,Qpu) 

+ AE{Qid,Qpd) (18) 



The next step is to select out the points that satisfy the 
constraint ( 16 ) and find the values of the pocket sizes for 



which the energy is minimized. 

For the undoped case, it turns out that the energy is 
minimized when the pockets in the linear band are zero, 
while the pockets in the parabolic band have a nonzero 
value. This is the result that one obtains if a monolayer 
and a bilayer are superimposed on each other. There is 
a priori no reason for this to be the case because in the 
exchange integrals there appear terms that are mixed in 
linear and parabolic pocket parameters. However, their 
contribution is too small to shift the equilibrium value 
of the pockets in the linear bands away from zero. In 
Fig. |4j we have plotted AE as function of Qi and Qp, 
where Qpu — —Qpd = Qp due to particle number con- 
servation and we have chosen Qi^ = —Qid = Qi- Since 
the spin of the electrons has no preferred direction, one 
sees two minima in Fig. |4] for Qi = and some fixed 
value of Qp = iQmin. The energy increases if the linear 
pocket is chosen to be different from zero, while tuning 
the parabolic pocket away from zero lowers the energy. 
Although AE is small (order of 1 meV per square mi- 
crometer), the equilibrium sizes of the pockets are sig- 
nificant (see Fig. [5| . The effect is comparable in magni- 
tude with the graphene bilayer. In Fig. |6]we display the 
equilibrium value for Qp as a function of the interaction 
strength g (for suspended graphene, g is estimated to be 
g « 2.3). The equilibrium value for the linear pocket 
sizes is zero for this range of the interaction strength. 
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The doped case in trilayer graphene is more subtle than 
in either a monolayer or a bilayer. For monolayer and bi- 
layers one can dope the system (with electrons or holes) 
and the bands (spin up and spin down) will fill up to some 
Fermi energy, corresponding with this particular doping 
level. This will be the noninteracting ground state for 
the doped system. In trilayer graphene this is not the 
case. If one dopes a graphene trilayer such that both the 
linear and the parabolic band are filled up to some Fermi 
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Figure 6. (Color online) Qmin, which is the equilibrium value 
of Qp as a function of the dimensionless interaction strength 
<?■ Qmin is measured in units of A. The line is a polynomial 
fit to eighth order in g. 



energy Ep, it turns out that due to kinetic energy consid- 
erations, this is not a stable state. The kinetic energy is 
minimized when the parabolic band is filled up differently 
than the linear band. Alternatively, since for a physical 
system the Fermi energy has a well defined value, one can 
interpret this result as a shift of the linear and parabolic 
energy bands with respect to each other. For our discus- 
sion, it is more natural to keep the intersection points of 
the bands in place and, as a consequence, use different 
Fermi energies for the parabolic and linear bands. By 
choosing this interpretation, we allow ourselves to use 
the formalism developed in the previous section. 

Since the kinetic energy cost of filling up the linear 
band goes as « k^, this costs more energy than filling 
up the parabolic band, for which the energy cost goes as 
w (recall that we work in dimensionless units, such 
that k < 1). Let us define kp^'"^'^ as the momentum 
to which the linear/parabolic spin up/down band fills up 
when the kinetic energy is minimized. When there is 
no interaction present the bands will be spin degenerate. 
Furthermore, we can use the same formalism as for the 
undoped case. The difference is that, for g = 0, the 
pocketsizes of the bands are equal to Q^^^ = kp^'"^'^. 



Hence, the constraint ( 16 ) now reads 



Slu- 



9L 

An 



Sld 



Mr 



9li + 

An 
P 2n - 



Qpu 

An 



^pd ■ 



Q^ 



pd 

An 



(19) 



where n is the doping level and s^^^ is the sign of Q^/p- 
To determine the values of one can vary the filling 
of the bands respecting the constraint and determine for 
which configuration the kinetic energy is minimized. One 
can show that Q° << Qp (see Fig. m). In fact, the res- 
olution we use for calculating the integrals is such that 
Q^' = 0. Note that, although Q^* << Q^, the single par- 
ticle energies associated with these momenta are of the 
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Figure 7. (Color online) Plot of Q1 (blue/black dots) and Qp 
(grey/grey dots) in units of A as a function of doping in units 
of hr^ . The red and blue lines mark the interval in which we 
have chosen our datapoints. The yellow line is a plot of Qp 
assuming that Qj' = 0. 



same order of magnitude. The linear band is filled to 
higher energies than the parabolic one, since the latter is 
very flat. However, in our discussion this is not relevant 
because the energies we calculate depend only on mo- 
menta and the fact that = in our formalism barely 
changes the results. Furthermore, if the effect of interac- 
tions on the linear pockets would be such that it would 
make them larger than the threshold value in Fig. [TJ we 
would be able to detect it. In the language we proposed 
in the introduction, this would be a band ferromagnetic 
state as the bands filled up to different energies, but have 
no net magnetization. 

In the doped case, the reference state with respect 
to which we compute energy differences has nonzero ki- 
netic and exchange energies, E^:^^ = E]^m{Q^iQ^) and 
E'^y. = 2Ecx{Q^, Qp)- One is now ready to vary the pocket 
parameters, compute the energies, apply the constraint 
(19), and find the configuration that minimizes the en- 



ergy 

A-E = i?kin(Q;n, Qpn) + E]^in{Qldi Qpd) ~ -^'kin 

+ Ec^{Qi 

Qpu) ~i~ E^^i^Qifij Qpd) ^cx' 

The result will depend on the value of the interaction pa- 
rameter g. If the graphene trilayer is doped, the system 
can still relax into a ferromagnetic state, but a critical 
interaction strength is needed. This critical value of the 
interaction increases with doping, as it can be seen in 
Fig. [sj The linear bands stay empty (up to our resolu- 
tion) and the parabolic pockets exhibit a discontinuous 
jump, indicating a first-order phase transition. This state 
is both band ferromagnetic, as well as spin ferromagnetic. 
Note that the jump is such that in one of the parabolic 
bands hole pockets will occur. 

So far we have looked only at configurations in which 
the pocket sizes are small. Although for the doped case 
the phase transition is first order, the pocket sizes are 
small and it is known that in monolayer graphene another 



first order transition occurs as the interaction strength 
exceeds some critical value {gc ~ 5.3 for undoped mono- 
layer graphene) ^ This transition is to a phase in which 
the monolayer has maximal magnetization. Since, for 
some purposes, one can regard a trilayer as a combina- 
tion of monolayer and bilayer graphene, it is natural to 
look for this transition in a graphene trilayer. Although 
this transition is theoretically present, we conclude that 
it can not been seen in any realistic experiment, because 
the critical coupling is out of any experimental range 

{gc > 200). 



IV. CONCLUSIONS 

In this paper, we have determined the ground state of 
trilayer graphene accounting for the long range Coulomb 
interaction. We used a formalism in which we could 
treat electron- and hole-like pockets on the same foot- 
ing. This allowed us to vary the four pocket parameters 
(linear/parabolic and spin up/down) to obtain a large 
dataset. We have chosen the discrete points to lie on a 
square root profile, so that we had many points that sat- 



isfied the constraint ( 16 1 for the undoped system or ( 19 ) 
for the doped one. 

For the undoped trilayer, we found that the energy is 
minimized for a configuration in which the linear bands 
are empty and an electron and a hole pocket occur in 
the spin up and spin down parabolic bands, which is 
a spin-ferromagnetic state (Fig. [2^). Since there is no 
preferred direction for the spin, this state is doubly de- 
generate (Fig. |4]) . The pockets increase in size when the 
interaction is tuned to higher values. They are only zero 
when the interaction vanishes, see Fig. [6j 

The doped trilayer is more subtle, since the noninter- 
acting case is already a band-ferromagnetic state in which 



n(10"') 




Figure 8. (Color online) Phase diagram, doping (n) versus 
interaction strength (g). The doping is dimensionless but can 
be converted to experimental units (cm~^) through multipli- 
cation with A'^. There is a first order phase transition from the 
ferromagnetic state (FM) to the normal state (N) as doping 
is increased. 
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the bands (linear/parabolic) fill up differently (Fig. [2]3). 
We named it a "band-ferromagnetic" state due to the fi- 
nite polarization in the pseudo-spin degree of freedom as- 
sociated with parabolic/linear bands. Although in phys- 
ical systems the bands will shift with respect to each 
other, resulting in a well defined Fermi energy, we chose 
to keep the bands fixed and let the bands fill up dif- 
ferently. This gave us two Fermi momenta (kp^) and 

Fermi energies {e'-J.''). Although E^p < E'-p, the parabohc 
band is much flatter than the linear one and therefore 
kp » kp. Our resolution was such that kp = 0, but 
this simplification will not affect the results. If the linear 
pockets exceed the threshold value given by the blue line 
in Fig. [7] for some value of the interaction strength we 
would have detected this. It turned out, however, that 
the linear bands stay empty for all doping levels that we 
considered. Furthermore, we saw a transition to a spin- 
ferromagnetic state. In contrast with the undoped case, 
this state is the ground state only if the coupling exceeds 
some critical value, which on its turn increases with dop- 
ing. The doping versus interaction strength phase dia- 
gram is shown in Fig. [8j The phase transition from the 
normal state (N) to a magnetic state (FM) is first order, 
i.e. the pocket size jumps discontinuously and the mag- 
netization also exhibits a jump to some nonzero value. 
Note that this magnetic state is both spin-ferromagnetic 
and band-ferromagnetic, since the bands fill up to differ- 
ent energies. 

We have also looked for a phase transition to a 
maximally magnetized state, as observed in monolayer 
graphene. We do not find such a transition for any inter- 



action strength that would be experimentally achievable. 

Although the graphene trilayer exhibits some features 
of both monolayer and bilayer graphene, it is an interest- 
ing system on itself and more complex than either of the 
two. The interplay between the filling of the linear and 
parabolic bands gives rise to many more possible config- 
urations of the pocket parameters. For example, already 
in the noninteracting groundstate of the doped trilayer 
the bands are shifted with respect to each other. 

It would be interesting to measure this spectrum in 
experiments using, for example, angle resolved photo- 
emission spectroscopy (ARPES). Long range Coulomb 
interactions can give rise to a ferromagnetic groundstate 
as it does in bilayer graphene, but will not affect the 
linear bands. The first order transition as seen in mono- 
layer graphene is not present as a result of interactions 
between the different bands. 

We are aware that next-nearest-neighbor hopping pa- 
rameters have effects on the energy spectrum that are of 
comparable magnitude as the effect we describe here.l^ 
However, if the system is sufficiently doped this will not 
alter our results. For the undoped case the results may 
be slightly altered, but our results could definitely be 
used as a starting point to investigate the full parameter 
model in more detail. 
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